#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import opengate as gate
from scipy.spatial.transform import Rotation
from pathlib import Path

# this first line is required at the beginning of all scripts
if __name__ == "__main__":
    sim = gate.Simulation()

    sim.g4_verbose = False
    sim.g4_verbose_level = 1
    sim.visu = False
    sim.visu_type = "vrml"
    sim.random_engine = "MersenneTwister"
    sim.random_seed = "auto"
    sim.number_of_threads = 4
    sim.progress_bar = True
    sim.output_dir = "./output"
    data_path = Path("data")

    # units
    m = gate.g4_units.m
    cm = gate.g4_units.cm
    nm = gate.g4_units.nm
    cm3 = gate.g4_units.cm3
    keV = gate.g4_units.keV
    MeV = gate.g4_units.MeV
    mm = gate.g4_units.mm
    Bq = gate.g4_units.Bq
    gcm3 = gate.g4_units.g / cm3

    # world
    world = sim.world
    world.size = [3 * m, 3 * m, 3 * m]
    world.material = "G4_AIR"

    # ct image
    patient = sim.add_volume("Image", "patient")
    patient.image = data_path / "patient-2mm.mhd"
    if sim.visu:
        # if the visualisation is enabled we load a very crude ct
        # otherwise the visualisation is too slow
        patient.image = data_path / "patient-20mm.mhd"
    patient.material = "G4_AIR"  # material used by default
    f1 = data_path / "Schneider2000MaterialsTable.txt"
    f2 = data_path / "Schneider2000DensitiesTable.txt"
    tol = 0.2 * gcm3
    (
        patient.voxel_materials,
        materials,
    ) = gate.geometry.materials.HounsfieldUnit_to_material(sim, tol, f1, f2)
    patient.rotation = Rotation.from_euler("x", 90, degrees=True).as_matrix()

    # CBCT gantry source
    gantry = sim.add_volume("Box", "CBCT_gantry")
    gantry.size = [0.2 * m, 0.2 * m, 0.2 * m]
    gantry.material = "G4_AIR"
    gantry.color = [0, 1, 1, 1]
    gantry.translation = [1060 * mm, 0, 0]

    # CBCT detector plane
    detector_plane = sim.add_volume("Box", "CBCT_detector_plane")
    detector_plane.size = [10 * mm, 409.6 * mm, 409.6 * mm]
    detector_plane.material = "G4_AIR"
    detector_plane.color = [1, 0, 0, 1]
    detector_plane.translation = [-536 * mm, 0, 0]

    # actor
    detector_actor = sim.add_actor("FluenceActor", "detector_actor")
    detector_actor.attached_to = detector_plane
    detector_actor.output_filename = "cbct.mhd"
    detector_actor.spacing = [10 * mm, 5 * mm, 5 * mm]
    detector_actor.size = [1, 100, 100]

    # physics
    sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option1"
    sim.physics_manager.set_production_cut("world", "all", 10 * mm)

    # source
    source = sim.add_source("GenericSource", "mysource")
    source.attached_to = gantry.name
    source.particle = "gamma"
    source.energy.mono = 60 * keV
    source.position.type = "box"
    source.position.size = [1 * nm, 16 * mm, 16 * mm]
    source.direction.type = "focused"
    source.direction.focus_point = [gantry.translation[0] - 60 * mm, 0, 0]
    source.n = 500000 / sim.number_of_threads

    # statistics
    stats = sim.add_actor("SimulationStatisticsActor", "stats")
    stats.track_types_flag = True
    stats.output_filename = "stats.txt"
    stats.write_to_disk = True

    # run
    sim.run()

    # print output statistics
    print(stats)
